load packages

library(tidyverse)
library(knitr)

define palettes

pal_control = c("#2E86B4", "#F2B827")
pal_outcome = wesanderson::wes_palette("Zissou1", 3, "continuous")
pal_condition = wesanderson::wes_palette("Zissou1", 6, "continuous")
pal_condition8 = c("#3B9AB2", "#357382", "#E86900", "#7FB8BA", "#ED4100", "#CC1400", "#EECC20", "#E1B002")
pal_condition8 = c("#CC1400", "#ED4100", "#E86900", "#357382", "#3B9AB2", "#7FB8BA", "#E1B002", "#EECC20")

source data

Note that % fat variable was incorrect in 2012_FDES –> replaced it with the data from the EMA dataframe

source("load_data.R")
saveRDS(ind_diffs, "ind_diffs.RDS")
saveRDS(ema, "ema.RDS")

ind_diffs = ind_diffs %>%
  mutate(bmi_z = as.numeric(scale(bmi)),
         fat_z = as.numeric(scale(fat)))

standardize brain data

Windsorizing observations +/- 3 SD from the mean to 3 –> otherwise we lose these subs in the SCA analyses
Z-score within ROI and session for ROIs
Z-score within map, test (association or uniformity), mask (masked or unmasked) and session for PEVs

# roi distributions
betas %>%
  filter(session == "all") %>%
  select(subjectID, con, condition, control, roi, process, meanPE) %>%
  unique() %>%
  filter(!is.na(roi)) %>%
  ggplot(aes(meanPE, fill = roi)) +
  geom_density(color = NA) +
  facet_wrap(~roi, ncol = 4) +
  theme_minimal(base_size = 14) +
  theme(legend.position = "none")

# dot product distribution
dots %>%
  filter(session == "all") %>%
  select(subjectID, con, condition, control, process, test, dotProduct) %>%
  unique() %>%
  ggplot(aes(dotProduct, fill = test)) +
  geom_density(color = NA) +
  facet_wrap(~process, ncol = 3, scales = "free") +
  theme_minimal(base_size = 14) +
  theme(legend.position = c(.85, .2))

dots %>%
  filter(session == "all") %>%
  select(subjectID, con, condition, control, process, test, dotProduct, mask) %>%
  unique() %>%
  ggplot(aes(dotProduct, fill = mask, alpha = .5)) +
  geom_density(color = NA) +
  facet_wrap(~process, ncol = 3, scales = "free") +
  theme_minimal(base_size = 14) +
  theme(legend.position = c(.85, .2))

# standardize and winsorize
betas_std = betas %>%
  group_by(roi, session) %>%
  mutate(meanPE_std = scale(meanPE, center = TRUE, scale = TRUE),
         outlier = ifelse(meanPE_std > 3, "yes",
                   ifelse(meanPE_std < -3, "yes", "no")),
         meanPE_std = ifelse(meanPE_std > 3, 3,
                      ifelse(meanPE_std < -3, -3, meanPE_std))) %>%
  ungroup()

dots_std = dots %>%
  group_by(map, test, mask, session) %>%
  mutate(dotProduct_std = scale(dotProduct, center = TRUE, scale = TRUE),
         outlier = ifelse(dotProduct_std > 3, "yes",
                   ifelse(dotProduct_std < -3, "yes", "no")),
         dotProduct_std = ifelse(dotProduct_std > 3, 3,
                          ifelse(dotProduct_std < -3, -3, dotProduct_std))) %>%
  ungroup()

summarize outliers

# summarize
betas_std %>%
  filter(session == "all") %>%
  group_by(outlier) %>%
  summarize(n = n()) %>%
  spread(outlier, n) %>%
  mutate(percent = round((yes / (yes + no)) * 100, 2))

dots_std %>%
  filter(session == "all") %>%
  group_by(outlier, test, mask) %>%
  summarize(n = n()) %>%
  spread(outlier, n) %>%
  mutate(percent = round((yes / (yes + no)) * 100, 2))

# remove outlier column
betas_std = betas_std %>%
  select(-outlier)

dots_std = dots_std %>%
  select(-outlier)

merge data

dataset = full_join(betas_std, dots_std, by = c("subjectID", "con", "process", "condition", "control", "session")) %>%
  filter(!condition == "social" & mask == "unmasked") %>%
  mutate(subjectID = as.character(subjectID)) %>%
  left_join(., ind_diffs, by = "subjectID") %>%
  ungroup() %>%
  mutate(subjectID = as.factor(subjectID),
         condition = as.factor(condition),
         control = as.factor(control),
         roi = as.factor(roi),
         process = as.factor(process),
         test = as.factor(test))

saveRDS(dataset, "full_dataset.RDS")
saveRDS(betas_std, "betas_dataset.RDS")
saveRDS(dots_std, "dots_dataset.RDS")

tidy data for plotting

Select cons generated across all sessions
Calculate mean across ROIs within each psychological process
Calculate balance score ROI = cognitive control - reward

ema_tidy = ema %>%
  select(subjectID, enact_prop) %>%
  unique() %>%
  filter(!is.na(subjectID)) %>%
  mutate(enact_prop_z = as.numeric(scale(enact_prop)))

dots_plot = dots_std %>%
  filter(!condition == "social" & mask == "unmasked" & session == "all") %>%
  select(-c(map, con, mask, session, dotProduct)) %>%
  left_join(., ind_diffs) %>%
  select(-c(sample, DBIC_ID, age, gender)) %>%
  left_join(., ema_tidy, by = "subjectID")

betas_plot = betas_std %>%
  filter(!condition == "social" & session == "all") %>%
  group_by(subjectID, process, condition, control) %>%
  mutate(meanProcessPEstd = mean(meanPE_std, na.rm = TRUE)) %>%
  select(-c(con, xyz, roi, meanPE, sdPE, meanPE_std)) %>%
  unique() %>%
  spread(process, meanProcessPEstd) %>%
  gather(process, meanProcessPEstd, cognitive_control, reward, value) %>%
  select(-c(session)) %>%
  left_join(., ind_diffs) %>%
  select(-c(sample, DBIC_ID, age, gender)) %>%
  ungroup() %>%
  left_join(., ema_tidy, by = "subjectID")

combined_plot = dots_plot %>%
  rename("value_std" = dotProduct_std) %>%
  mutate(test = ifelse(process == "craving_regulation", "neural signature", test)) %>%
  select(subjectID, process, test, condition, control, value_std, fat, bmi, enact_prop) %>%
  bind_rows(betas_plot %>%
              rename("value_std" = meanProcessPEstd) %>%
              mutate(test = "ROI") %>%
              select(-contains("_z")))

tidy data for modeling

model_data = dataset %>%
  filter(test == "association" & session == "all" & control == "nature" & condition == "food") %>%
  unique() %>%
  group_by(process, subjectID) %>%
  mutate(meanProcessPEstd = mean(meanPE_std, na.rm = TRUE),
         processPE = paste0(process, "_ROI"),
         processPEV = paste0(process, "_PEV")) %>%
  ungroup() %>%
  select(-c(xyz, meanPE, meanPE_std, sdPE, dotProduct, map, process)) %>%
  spread(processPEV, dotProduct_std) %>%
  spread(processPE, meanProcessPEstd) %>%
  group_by(subjectID) %>%
  fill(contains("PEV"), .direction = "down") %>%
  fill(contains("PEV"), .direction = "up") %>%
  fill(contains("ROI"), .direction = "down") %>%
  fill(contains("ROI"), .direction = "up") %>%
  select(-c(roi, craving_ROI, craving_regulation_ROI, test, mask, age, gender)) %>%
  unique()

model_data_bmi = model_data %>%
  select(-contains("fat"), -restraint) %>%
  na.omit()

model_data_fat = model_data %>%
  select(-contains("bmi"), -restraint) %>%
  na.omit()

model_data_ema = model_data %>%
  select(-contains("fat"), -contains("bmi"), -restraint) %>%
  right_join(., ema_tidy, by = "subjectID") %>%
  na.omit()

visualize

DVs and dietary restraint

# source raincloud plot
source("https://gist.githubusercontent.com/benmarwick/2a1bb0133ff568cbe28d/raw/fb53bd97121f7f9ce947837ef1a4c65a73bffb3f/geom_flat_violin.R")

# plot
ind_diffs %>%
  left_join(., ema_tidy, by = "subjectID") %>%
  select(subjectID, bmi, fat, enact_prop, restraint) %>%
  rename("BMI" = bmi,
         "body fat percentage" = fat,
         "enactment" = enact_prop,
         "dietary restraint" = restraint) %>%
  gather(var, val, -subjectID) %>%
  mutate(var = factor(var, levels = c("BMI", "body fat percentage", "enactment", "dietary restraint"))) %>%
  ggplot(aes("", val, fill = var, color = var)) +
    geom_flat_violin(position = position_nudge(x = .1, y = 0), color = FALSE) +
    geom_boxplot(width = .1, outlier.shape = NA, color = "black", fill = "white") +
    geom_point(position = position_jitter(width = .05), size = .5, alpha = .75) +
    facet_wrap(~var, scales = "free", nrow = 1) +
    scale_fill_manual(values = c(pal_outcome, pal_condition[3])) +
    scale_color_manual(values = c(pal_outcome, pal_condition[3])) +
    scale_x_discrete(expand = c(.1, 0)) +
    coord_flip() +
    labs(y = "\nvalue", x = "") +
    theme_minimal(base_size = 14) +
    theme(legend.position = "none")

combined

combined_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "nature")),
         process = gsub("_", " ", process),
         test = factor(test, levels = c("neural signature", "ROI", "association", "uniformity"))) %>%
  ggplot(aes(control, value_std, fill = condition)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "bar", position = position_dodge(.9)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", position = position_dodge(.9), width = 0) +
  facet_wrap(~ test + process, nrow = 3) +
  scale_fill_manual(name = "", values = pal_condition) +
  labs(x = "\ncontrol condition", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

combined_plot %>%
  filter(!condition == "nature") %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "nature")),
         process = gsub("_", " ", process),
         process = ifelse(!test == "ROI", sprintf("%s PEV", process), sprintf("%s ROI", process)),
         process = factor(process, levels = c("reward ROI", "reward PEV", "craving PEV",
                                              "cognitive control ROI", "cognitive control PEV", "craving regulation PEV",
                                              "value ROI", "value PEV")),
         test = ifelse(!test %in% c("association", "uniformity"), "none", test),
         test = factor(test, levels = c("none", "association", "uniformity"))) %>%
  ggplot(aes(control, value_std, fill = process)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "bar", position = position_dodge(.9)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", position = position_dodge(.9), width = 0) +
  facet_grid(test~condition) +
  scale_fill_manual(name = "", values = pal_condition8) +
  labs(x = "\ncontrol condition", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

ROIs

means

Mean parameter estimates as a function of process, stimulus type, and control condition

betas_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "nature"))) %>%
  ggplot(aes(condition, meanProcessPEstd, fill = control)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "bar", position = position_dodge(.9)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", position = position_dodge(.9), width = 0) +
  facet_grid(~process) +
  scale_fill_manual(name = "", values = pal_control) +
  labs(x = "\nstimulus", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

betas_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "nature"))) %>%
  ggplot(aes(control, meanProcessPEstd, fill = condition)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "bar", position = position_dodge(.9)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", position = position_dodge(.9), width = 0) +
  facet_grid(~process) +
  scale_fill_manual(name = "", values = pal_condition) +
  labs(x = "\ncontrol condition", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

correlations

Correlations between variables and BMI, % fat, or EMA enactment proportion

betas_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "nature"))) %>%
  ggplot(aes(bmi, meanProcessPEstd, color = condition)) +
  geom_point(alpha = .2) +
  geom_smooth(method = "lm", alpha = .25) +
  facet_grid(control~process) +
  scale_color_manual(name = "", values = pal_condition) +
  labs(x = "\nBMI", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

betas_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "nature"))) %>%
  ggplot(aes(fat, meanProcessPEstd, color = condition)) +
  geom_point(alpha = .2) +
  geom_smooth(method = "lm", alpha = .25) +
  facet_grid(control~process) +
  scale_color_manual(name = "", values = pal_condition) +
  labs(x = "\nbody composition (% fat)", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

betas_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "nature"))) %>%
  ggplot(aes(enact_prop, meanProcessPEstd, color = condition)) +
  geom_point(alpha = .2) +
  geom_smooth(method = "lm", alpha = .25) +
  facet_grid(control~process) +
  scale_color_manual(name = "", values = pal_condition) +
  labs(x = "\ndaily enactment proportion", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

PEVs

means

Mean expression as a function of process, stimulus type, control condition, and test

dots_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "nature"))) %>%
  ggplot(aes(condition, dotProduct_std, fill = control)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "bar", position = position_dodge(.9)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", position = position_dodge(.9), width = 0) +
  facet_grid(test~process) +
  scale_fill_manual(name = "", values = pal_control) +
  labs(x = "\nstimulus", y = "mean expression\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

dots_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "nature"))) %>%
  ggplot(aes(control, dotProduct_std, fill = condition)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "bar", position = position_dodge(.9)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", position = position_dodge(.9), width = 0) +
  facet_grid(test~process) +
  scale_fill_manual(name = "", values = pal_condition) +
  labs(x = "\ncontrol condition", y = "mean expression\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

correlations

Correlations between variables and BMI, % fat, or EMA enactment proportion

dots_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature"))) %>%
  ggplot(aes(bmi, dotProduct_std, color = condition)) +
  geom_point(alpha = .2) +
  geom_smooth(method = "lm", alpha = .25) +
  facet_grid(control + test ~ process) +
  scale_color_manual(name = "", values = pal_condition) +
  labs(x = "\nBMI", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

dots_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature"))) %>%
  ggplot(aes(fat, dotProduct_std, color = condition)) +
  geom_point(alpha = .2) +
  geom_smooth(method = "lm", alpha = .25) +
  facet_grid(control + test ~ process) +
  scale_color_manual(name = "", values = pal_condition) +
  labs(x = "\nbody composition (% fat)", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

dots_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature"))) %>%
  ggplot(aes(enact_prop, dotProduct_std, color = condition)) +
  geom_point(alpha = .2) +
  geom_smooth(method = "lm", alpha = .25) +
  facet_grid(control + test ~ process) +
  scale_color_manual(name = "", values = pal_condition) +
  labs(x = "\nbody composition (% fat)", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

correlations among ROIs and PEVs

PEVs v ROIs

dots_plot %>%
  unite(var, process, test, condition, control, sep = "_") %>%
  mutate(var = sprintf("dot_%s", var)) %>%
  spread(var, dotProduct_std) %>%
  left_join(., betas_plot) %>%
  unite(var, process, condition, control, sep = "_") %>%
  mutate(var = sprintf("roi_%s", var)) %>%
  spread(var, meanProcessPEstd) %>%
  select(-subjectID, -contains("bmi"), -contains("fat"), -contains("restraint"), -contains("enact_prop")) %>%
  cor(., use = "pairwise.complete.obs") %>%
  reshape2::melt() %>%
  filter(!grepl("dot", Var1) & !grepl("roi", Var2)) %>%
  ggplot(aes(Var1, Var2, fill = value)) +
    geom_tile(color = "white") +
    scale_fill_gradientn(name = "correlation\n", colors = c(pal_condition[1], "white", pal_condition[6]), limits = c(-1, 1), breaks = c(-1, 0, 1)) + 
    geom_text(aes(label = round(value, 1)), size = 2) +
    labs(x = "", y = "") + 
    theme(legend.text = element_text(size = 8),
          legend.position = "top",
          legend.title = element_blank(),
          axis.text = element_text(color = "black"),
          axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
          axis.line = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          panel.background = element_blank())

ROIs

dots_plot %>%
  unite(var, process, test, condition, control, sep = "_") %>%
  mutate(var = sprintf("dot_%s", var)) %>%
  spread(var, dotProduct_std) %>%
  left_join(., betas_plot) %>%
  unite(var, process, condition, control, sep = "_") %>%
  mutate(var = sprintf("roi_%s", var)) %>%
  spread(var, meanProcessPEstd) %>%
  select(-subjectID) %>%
  cor(., use = "pairwise.complete.obs") %>%
  reshape2::melt() %>%
  filter(grepl("roi", Var1) & grepl("roi", Var2)) %>%
  ggplot(aes(Var1, Var2, fill = value)) +
    geom_tile(color = "white") +
    scale_fill_gradientn(name = "correlation\n", colors = c(pal_condition[1], "white", pal_condition[6]), limits = c(-1, 1), breaks = c(-1, 0, 1)) + 
    geom_text(aes(label = round(value, 1)), size = 2) +
    labs(x = "", y = "") + 
    theme(legend.text = element_text(size = 8),
          legend.position = "top",
          legend.title = element_blank(),
          axis.text = element_text(color = "black"),
          axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
          axis.line = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          panel.background = element_blank())

PEVs

dots_plot %>%
  unite(var, process, test, condition, control, sep = "_") %>%
  mutate(var = sprintf("dot_%s", var)) %>%
  spread(var, dotProduct_std) %>%
  left_join(., betas_plot) %>%
  unite(var, process, condition, control, sep = "_") %>%
  mutate(var = sprintf("roi_%s", var)) %>%
  spread(var, meanProcessPEstd) %>%
  select(-subjectID) %>%
  cor(., use = "pairwise.complete.obs") %>%
  reshape2::melt() %>%
  filter(grepl("dot", Var1) & grepl("dot", Var2)) %>%
  ggplot(aes(Var1, Var2, fill = value)) +
    geom_tile(color = "white") +
    scale_fill_gradientn(name = "correlation\n", colors = c(pal_condition[1], "white", pal_condition[6]), limits = c(-1, 1), breaks = c(-1, 0, 1)) + 
    geom_text(aes(label = round(value, 1)), size = 2) +
    labs(x = "", y = "") + 
    theme(legend.text = element_text(size = 8),
          legend.position = "top",
          legend.title = element_blank(),
          axis.text = element_text(color = "black"),
          axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
          axis.line = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          panel.background = element_blank())

descriptives

means and sds

model_data %>%
  left_join(., model_data_ema) %>%
  gather(Variable, val, contains("PEV"), contains("ROI"), bmi, fat, enact_prop, restraint) %>%
  filter(!Variable == "balance_ROI") %>%
  group_by(Variable) %>%
  summarize(M = mean(val, na.rm = TRUE),
            SD = sd(val, na.rm = TRUE)) %>%
  filter(!is.na(Variable)) %>%
  mutate(Variable = gsub("_", " ", Variable),
         Variable = ifelse(Variable == "bmi", "BMI",
                    ifelse(Variable == "fat", "Body fat percentage",
                    ifelse(Variable == "enact prop", "Enactment",
                    ifelse(Variable == "restraint", "Dietary restraint", Hmisc::capitalize(Variable))))),
         Variable = factor(Variable, levels = c("BMI", "Body fat percentage", "Enactment",
                                                "Reward ROI", "Reward PEV", "Craving PEV",
                                                "Cognitive control ROI", "Cognitive control PEV", "Craving regulation PEV",
                                                "Value ROI", "Value PEV",
                                                "Dietary restraint"))) %>%
  arrange(Variable) %>%
  kable(format = "pandoc", digits = 2)

Variable	M	SD
BMI	23.19	3.36
Body fat percentage	25.43	8.00
Enactment	0.47	0.22
Reward ROI	0.06	0.50
Reward PEV	0.07	0.77
Craving PEV	-0.13	0.76
Cognitive control ROI	-0.13	0.50
Cognitive control PEV	-0.11	0.78
Craving regulation PEV	0.40	0.67
Value ROI	0.40	0.68
Value PEV	0.28	0.77
Dietary restraint	16.90	8.10

correlations

cor_mat = model_data %>%
  left_join(., model_data_ema) %>%
  select(subjectID, bmi, fat, enact_prop, contains("ROI"), contains("PEV"), restraint) %>%
  gather(var, val, contains("PEV"), contains("ROI"), bmi, fat, enact_prop, restraint) %>%
  mutate(var = gsub("_", " ", var),
         var = ifelse(var == "bmi", "BMI",
               ifelse(var == "fat", "Body fat percentage",
               ifelse(var == "enact prop", "Enactment",
               ifelse(var == "restraint", "Dietary restraint", Hmisc::capitalize(var))))),
         var = factor(var, levels = c("BMI", "Body fat percentage", "Enactment",
                                      "Reward ROI", "Reward PEV", "Craving PEV",
                                      "Cognitive control ROI", "Cognitive control PEV", "Craving regulation PEV",
                                      "Value ROI", "Value PEV",
                                      "Dietary restraint"))) %>%
  arrange(var) %>%
  spread(var, val) %>%
  ungroup() %>%
  select(-subjectID) %>%
  as.matrix()

cors = Hmisc::rcorr(cor_mat)

p = cors$P %>%
  as.data.frame() %>%
  rownames_to_column() %>%
  gather(var, p_value, -rowname) %>%
  mutate(code = ifelse(between(p_value, .05, .1), "†",
                ifelse(between(p_value, .01, .05), "*", 
                ifelse(between(p_value, .001, .01), "**",
                ifelse(p_value < .001, "***", "")))))

cors$r %>%
  as.data.frame() %>%
  rownames_to_column() %>%
  gather(var, val, -rowname) %>%
  left_join(., p) %>%
  mutate(val = gsub("0.(.*)", ".\\1", sprintf("%.2f%s", val, code)),
         val = gsub("1..NA", "--", val),
         rowname = factor(rowname, levels = c("BMI", "Body fat percentage", "Enactment",
                                      "Reward ROI", "Reward PEV", "Craving PEV",
                                      "Cognitive control ROI", "Cognitive control PEV", "Craving regulation PEV",
                                      "Value ROI", "Value PEV",
                                      "Dietary restraint"))) %>%
  select(-p_value, -code) %>%
  spread(var, val) %>%
  rename(" " = rowname) %>%
  select(" ", "BMI", "Body fat percentage", "Enactment",
         "Reward ROI", "Reward PEV", "Craving PEV",
         "Cognitive control ROI", "Cognitive control PEV", "Craving regulation PEV",
         "Value ROI", "Value PEV",
         "Dietary restraint") %>%
  kable(format = "pandoc")

	BMI	Body fat percentage	Enactment	Reward ROI	Reward PEV	Craving PEV	Cognitive control ROI	Cognitive control PEV	Craving regulation PEV	Value ROI	Value PEV	Dietary restraint
BMI	–	.66***	-.13	-.12†	-.02	-.02	-.02	-.06	.10	.06	.00	.24***
Body fat percentage	.66***	–	-.20†	-.07	-.07	-.07	-.03	-.06	.02	.01	-.04	.36***
Enactment	-.13	-.20†	–	.21*	.31**	.21*	.01	.05	.03	.23*	.29**	.14
Reward ROI	-.12†	-.07	.21*	–	.62***	.47***	.31***	.41***	.08	.30***	.60***	-.02
Reward PEV	-.02	-.07	.31**	.62***	–	.81***	.47***	.59***	.19**	.58***	.91***	.01
Craving PEV	-.02	-.07	.21*	.47***	.81***	–	.45***	.54***	.24***	.51***	.76***	.10
Cognitive control ROI	-.02	-.03	.01	.31***	.47***	.45***	–	.91***	.21***	.06	.30***	-.15*
Cognitive control PEV	-.06	-.06	.05	.41***	.59***	.54***	.91***	–	.23***	.16**	.40***	-.13*
Craving regulation PEV	.10	.02	.03	.08	.19**	.24***	.21***	.23***	–	.28***	.15*	.02
Value ROI	.06	.01	.23*	.30***	.58***	.51***	.06	.16**	.28***	–	.65***	.14*
Value PEV	.00	-.04	.29**	.60***	.91***	.76***	.30***	.40***	.15*	.65***	–	.07
Dietary restraint	.24***	.36***	.14	-.02	.01	.10	-.15*	-.13*	.02	.14*	.07	–

t-tests

model_data %>%
  select(-contains("z"), -contains("age"), -contains("restraint"), -bmi, -fat) %>%
  ungroup() %>%
  select_if(is.numeric) %>%
  map_df(~ broom::tidy(t.test(.x, mu = 0, na.action = "na.omit")), .id = "x") %>%
  rename("Variable" = x,
         "t" = statistic,
         "df" = parameter,
         "p" = p.value) %>%
  mutate(`Mdiff [95% CI]` = sprintf("%.02f [%.02f, %.02f]", estimate, conf.low, conf.high),
         p = ifelse(p < .001, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p))),
         Variable = gsub("_", " ", Variable),
         Variable = Hmisc::capitalize(Variable),
         Variable = factor(Variable, levels = c("Reward ROI", "Reward PEV", "Craving PEV",
                                                "Cognitive control ROI", "Cognitive control PEV", "Craving regulation PEV",
                                                "Value ROI", "Value PEV"))) %>%
  select(Variable, `Mdiff [95% CI]`, t, df, p) %>%
  arrange(Variable) %>%
  kable(format = "pandoc", digits = 2)

Variable	Mdiff [95% CI]	t	df	p
Reward ROI	0.06 [-0.00, 0.12]	1.90	261	.059
Reward PEV	0.07 [-0.02, 0.17]	1.51	261	.132
Craving PEV	-0.13 [-0.22, -0.04]	-2.74	261	.007
Cognitive control ROI	-0.13 [-0.19, -0.07]	-4.14	261	< .001
Cognitive control PEV	-0.11 [-0.20, -0.01]	-2.20	261	.029
Craving regulation PEV	0.40 [0.31, 0.48]	9.52	261	< .001
Value ROI	0.40 [0.31, 0.48]	9.38	261	< .001
Value PEV	0.28 [0.18, 0.37]	5.84	261	< .001

run cross-validated models

specify model variables

options(na.action = "na.fail")
data_ctrl = caret::trainControl(method = "repeatedcv", number = 5, repeats = 5)
set.seed(123)

BMI

ROI models

bmi_roi_full = lm(bmi_z ~ value_ROI + reward_ROI + cognitive_control_ROI, data = model_data_bmi)
(bmi_roi_mods = MuMIn::dredge(bmi_roi_full))

MuMIn::get.models(bmi_roi_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()

best_bmi_roi = caret::train(bmi_z ~ reward_ROI + value_ROI,
                     data = model_data_bmi,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

PEV models

bmi_pev_full = lm(bmi_z ~ value_PEV + reward_PEV + cognitive_control_PEV + craving_PEV + craving_regulation_PEV, data = model_data_bmi)
(bmi_pev_mods = MuMIn::dredge(bmi_pev_full))

MuMIn::get.models(bmi_pev_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()

best_bmi_pev = caret::train(bmi_z ~ craving_regulation_PEV,
                     data = model_data_bmi,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

combined

bmi_null = lm(bmi_z ~ 1, data = model_data_bmi)
summary(bmi_null)

## 
## Call:
## lm(formula = bmi_z ~ 1, data = model_data_bmi)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2100 -0.6542 -0.0558  0.4827  4.6115 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.01306    0.06326   0.206    0.837
## 
## Residual standard error: 1.006 on 252 degrees of freedom

bmi_combined = caret::train(bmi_z ~ reward_ROI + value_ROI + craving_regulation_PEV,
                     data = model_data_bmi,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

(bmi_aic = AIC(bmi_null, best_bmi_roi$finalModel, best_bmi_pev$finalModel, bmi_combined$finalModel) %>%
    rownames_to_column() %>%
    extract(rowname, "model", ".*bmi_(.*)\\$.*") %>%
    mutate(model = toupper(model),
           model = ifelse(is.na(model), "Null", model)) %>%
  arrange(AIC))

summarize the overall best fitting

# ROI
summary(best_bmi_roi$finalModel)

## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1355 -0.6606 -0.0899  0.4662  4.5052 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.03113    0.07236  -0.430   0.6674  
## reward_ROI  -0.30523    0.12990  -2.350   0.0196 *
## value_ROI    0.15583    0.09654   1.614   0.1077  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9976 on 250 degrees of freedom
## Multiple R-squared:  0.02502,    Adjusted R-squared:  0.01722 
## F-statistic: 3.207 on 2 and 250 DF,  p-value: 0.04214

best_bmi_roi

## Linear Regression 
## 
## 253 samples
##   2 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 5 times) 
## Summary of sample sizes: 202, 203, 202, 201, 204, 203, ... 
## Resampling results:
## 
##   RMSE       Rsquared    MAE      
##   0.9991722  0.03203531  0.7489974
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

sd(best_bmi_roi$resample$Rsquared)

## [1] 0.04890756

sd(best_bmi_roi$resample$RMSE)

## [1] 0.1078594

# combined
summary(bmi_combined$finalModel)

## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1994 -0.6654 -0.0886  0.5098  4.4311 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)  
## (Intercept)            -0.07322    0.07768  -0.942   0.3469  
## reward_ROI             -0.30708    0.12961  -2.369   0.0186 *
## value_ROI               0.11695    0.09989   1.171   0.2428  
## craving_regulation_PEV  0.14370    0.09792   1.467   0.1435  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9953 on 249 degrees of freedom
## Multiple R-squared:  0.03338,    Adjusted R-squared:  0.02173 
## F-statistic: 2.866 on 3 and 249 DF,  p-value: 0.03725

bmi_combined

## Linear Regression 
## 
## 253 samples
##   3 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 5 times) 
## Summary of sample sizes: 203, 201, 203, 202, 203, 203, ... 
## Resampling results:
## 
##   RMSE      Rsquared    MAE      
##   1.004276  0.04192202  0.7616355
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

sd(bmi_combined$resample$Rsquared)

## [1] 0.045541

sd(bmi_combined$resample$RMSE)

## [1] 0.1267501

% fat

ROI models

fat_roi_full = lm(fat_z ~ value_ROI + reward_ROI + cognitive_control_ROI, data = model_data_fat)
(fat_roi_mods = MuMIn::dredge(fat_roi_full))

MuMIn::get.models(fat_roi_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()

best_fat_roi = caret::train(fat_z ~ reward_ROI,
                     data = model_data_fat,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

PEV models

fat_pev_full = lm(fat_z ~ value_PEV + reward_PEV + cognitive_control_PEV + craving_PEV + craving_regulation_PEV, data = model_data_fat)
(fat_pev_mods = MuMIn::dredge(fat_pev_full))

MuMIn::get.models(fat_pev_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()

best_fat_pev = caret::train(fat_z ~ craving_PEV,
                     data = model_data_fat,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

combined

fat_null = lm(fat_z ~ 1, data = model_data_fat)

fat_combined = caret::train(fat_z ~ reward_ROI + craving_PEV,
                     data = model_data_fat,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

(fat_aic = AIC(fat_null, best_fat_roi$finalModel, best_fat_pev$finalModel, fat_combined$finalModel) %>%
    rownames_to_column() %>%
    extract(rowname, "model", ".*fat_(.*)\\$.*") %>%
    mutate(model = toupper(model),
           model = ifelse(is.na(model), "Null", model)) %>%
    arrange(AIC))

summarize the overall best fitting

summary(fat_null)

## 
## Call:
## lm(formula = fat_z ~ 1, data = model_data_fat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9531 -0.5945  0.1396  0.6457  2.4201 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.01679    0.06321   0.266    0.791
## 
## Residual standard error: 0.9994 on 249 degrees of freedom

EMA enactment proportion

ROI models

enact_prop_roi_full = lm(enact_prop_z ~ value_ROI + reward_ROI + cognitive_control_ROI, data = model_data_ema)
(enact_prop_roi_mods = MuMIn::dredge(enact_prop_roi_full))

MuMIn::get.models(enact_prop_roi_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()

best_enact_prop_roi = caret::train(enact_prop_z ~ value_ROI,
                     data = model_data_ema,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

PEV models

enact_prop_pev_full = lm(enact_prop_z ~ value_PEV + reward_PEV + cognitive_control_PEV + craving_PEV + craving_regulation_PEV, data = model_data_ema)
(enact_prop_pev_mods = MuMIn::dredge(enact_prop_pev_full))

MuMIn::get.models(enact_prop_pev_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()

best_enact_prop_pev = caret::train(enact_prop_z ~ reward_PEV,
                     data = model_data_ema,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

combined

enact_prop_null = lm(enact_prop_z ~ 1, data = model_data_ema)

enact_prop_combined = caret::train(enact_prop_z ~ value_ROI + reward_PEV,
                     data = model_data_ema,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

(enact_aic = AIC(enact_prop_null, best_enact_prop_roi$finalModel, best_enact_prop_pev$finalModel, enact_prop_combined$finalModel) %>%
    rownames_to_column() %>%
    extract(rowname, "model", ".*prop_(.*)\\$.*") %>%
    mutate(model = toupper(model),
           model = ifelse(is.na(model), "Null", model)) %>%
  arrange(AIC))

summarize the overall best fitting

summary(best_enact_prop_pev$finalModel)

## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.14970 -0.57410 -0.01054  0.53196  2.37344 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -0.06291    0.09858  -0.638  0.52488   
## reward_PEV   0.36683    0.11493   3.192  0.00192 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9573 on 95 degrees of freedom
## Multiple R-squared:  0.09686,    Adjusted R-squared:  0.08735 
## F-statistic: 10.19 on 1 and 95 DF,  p-value: 0.001917

best_enact_prop_pev

## Linear Regression 
## 
## 97 samples
##  1 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 5 times) 
## Summary of sample sizes: 78, 77, 78, 77, 78, 78, ... 
## Resampling results:
## 
##   RMSE       Rsquared   MAE      
##   0.9565301  0.1467801  0.7619614
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

sd(best_enact_prop_pev$resample$Rsquared)

## [1] 0.1304739

sd(best_enact_prop_pev$resample$RMSE)

## [1] 0.1571684

tables

bmi

bmi_fit = data.frame(model = c("ROI", "PEV", "COMBINED"),
                     r2 = c(mean(best_bmi_roi$resample$Rsquared), 
                            mean(best_bmi_pev$resample$Rsquared), 
                            mean(bmi_combined$resample$Rsquared)),
                     r2_sd = c(sd(best_bmi_roi$resample$Rsquared), 
                               sd(best_bmi_pev$resample$Rsquared), 
                               sd(bmi_combined$resample$Rsquared)),
                     RMSE = c(mean(best_bmi_roi$resample$RMSE), 
                              mean(best_bmi_pev$resample$RMSE), 
                              mean(bmi_combined$resample$RMSE)),
                     RMSE_sd = c(sd(best_bmi_roi$resample$RMSE), 
                                 sd(best_bmi_pev$resample$RMSE), 
                                 sd(bmi_combined$resample$RMSE)))

bmi_table = broom::tidy(best_bmi_roi$finalModel, conf.int = TRUE) %>% 
  mutate(model = "ROI", 
         mod_sig = ifelse(pf(summary(best_bmi_roi$finalModel)$fstatistic[1], 
                             summary(best_bmi_roi$finalModel)$fstatistic[2], 
                             summary(best_bmi_roi$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", "")) %>%
  bind_rows(broom::tidy(best_bmi_pev$finalModel, conf.int = TRUE) %>% 
              mutate(model = "PEV",
                     mod_sig = ifelse(pf(summary(best_bmi_pev$finalModel)$fstatistic[1], 
                             summary(best_bmi_pev$finalModel)$fstatistic[2], 
                             summary(best_bmi_pev$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", ""))) %>%
  bind_rows(broom::tidy(bmi_combined$finalModel, conf.int = TRUE) %>% 
              mutate(model = "COMBINED",
                     mod_sig = ifelse(pf(summary(bmi_combined$finalModel)$fstatistic[1], 
                             summary(bmi_combined$finalModel)$fstatistic[2], 
                             summary(bmi_combined$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", ""))) %>%
  bind_rows(broom::tidy(bmi_null, conf.int = TRUE) %>% mutate(model = "Null", mod_sig = "")) %>%
  rename("SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%.2f, %.2f]", estimate, conf.low, conf.high),
         p = ifelse(p < .001, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p))),
         term = gsub("\\(Intercept\\)", "Intercept", term),
         term = gsub("_", " ", term),
         term = Hmisc::capitalize(term)) %>%
  left_join(., bmi_aic) %>%
  left_join(., bmi_fit) %>%
  mutate(`r2 (SD)` = ifelse(!model == "Null", sprintf("%.2f (%.2f)", r2, r2_sd), "--"),
         `RMSE (SD)` = ifelse(!model == "Null", sprintf("%.2f (%.2f)", RMSE, RMSE_sd), "--")) %>%
  unite(model, model, mod_sig, sep = "") %>%
  select(model, AIC, `r2 (SD)`, `RMSE (SD)`, term, `b [95% CI]`, SE, t, p) %>%
  arrange(AIC) %>%
  mutate(model = gsub("COMBINED", "Combined", model),
         outcome = "BMI") %>%
    select(outcome, model, AIC, `r2 (SD)`, `RMSE (SD)`, term, `b [95% CI]`, SE, t, p)

bmi_table %>%
  kable(format = "pandoc", digits = 2)

outcome	model	AIC	r2 (SD)	RMSE (SD)	term	b [95% CI]	SE	t	p
BMI	Combined*	721.56	0.04 (0.05)	1.00 (0.13)	Intercept	-0.07 [-0.23, 0.08]	0.08	-0.94	.347
BMI	Combined*	721.56	0.04 (0.05)	1.00 (0.13)	Reward ROI	-0.31 [-0.56, -0.05]	0.13	-2.37	.019
BMI	Combined*	721.56	0.04 (0.05)	1.00 (0.13)	Value ROI	0.12 [-0.08, 0.31]	0.10	1.17	.243
BMI	Combined*	721.56	0.04 (0.05)	1.00 (0.13)	Craving regulation PEV	0.14 [-0.05, 0.34]	0.10	1.47	.144
BMI	ROI*	721.74	0.03 (0.05)	1.00 (0.11)	Intercept	-0.03 [-0.17, 0.11]	0.07	-0.43	.667
BMI	ROI*	721.74	0.03 (0.05)	1.00 (0.11)	Reward ROI	-0.31 [-0.56, -0.05]	0.13	-2.35	.020
BMI	ROI*	721.74	0.03 (0.05)	1.00 (0.11)	Value ROI	0.16 [-0.03, 0.35]	0.10	1.61	.108
BMI	PEV	723.45	0.02 (0.03)	0.99 (0.14)	Intercept	-0.05 [-0.19, 0.10]	0.07	-0.67	.507
BMI	PEV	723.45	0.02 (0.03)	0.99 (0.14)	Craving regulation PEV	0.16 [-0.03, 0.34]	0.09	1.64	.102
BMI	Null	724.15	–	–	Intercept	0.01 [-0.11, 0.14]	0.06	0.21	.837

% fat

fat_fit = data.frame(model = c("ROI", "PEV", "COMBINED"),
                     r2 = c(mean(best_fat_roi$resample$Rsquared), 
                            mean(best_fat_pev$resample$Rsquared), 
                            mean(fat_combined$resample$Rsquared)),
                     r2_sd = c(sd(best_fat_roi$resample$Rsquared), 
                               sd(best_fat_pev$resample$Rsquared), 
                               sd(fat_combined$resample$Rsquared)),
                     RMSE = c(mean(best_fat_roi$resample$RMSE), 
                              mean(best_fat_pev$resample$RMSE), 
                              mean(fat_combined$resample$RMSE)),
                     RMSE_sd = c(sd(best_fat_roi$resample$RMSE), 
                                 sd(best_fat_pev$resample$RMSE), 
                                 sd(fat_combined$resample$RMSE)))

fat_table = broom::tidy(best_fat_roi$finalModel, conf.int = TRUE) %>% 
  mutate(model = "ROI", 
         mod_sig = ifelse(pf(summary(best_fat_roi$finalModel)$fstatistic[1], 
                             summary(best_fat_roi$finalModel)$fstatistic[2], 
                             summary(best_fat_roi$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", "")) %>%
  bind_rows(broom::tidy(best_fat_pev$finalModel, conf.int = TRUE) %>% 
              mutate(model = "PEV",
                     mod_sig = ifelse(pf(summary(best_fat_pev$finalModel)$fstatistic[1], 
                             summary(best_fat_pev$finalModel)$fstatistic[2], 
                             summary(best_fat_pev$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", ""))) %>%
  bind_rows(broom::tidy(fat_combined$finalModel, conf.int = TRUE) %>% 
              mutate(model = "COMBINED",
                     mod_sig = ifelse(pf(summary(fat_combined$finalModel)$fstatistic[1], 
                             summary(fat_combined$finalModel)$fstatistic[2], 
                             summary(fat_combined$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", ""))) %>%
  bind_rows(broom::tidy(fat_null, conf.int = TRUE) %>% mutate(model = "Null", mod_sig = "")) %>%
  rename("SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%.2f, %.2f]", estimate, conf.low, conf.high),
         p = ifelse(p < .001, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p))),
         term = gsub("\\(Intercept\\)", "Intercept", term),
         term = gsub("_", " ", term),
         term = Hmisc::capitalize(term)) %>%
  left_join(., fat_aic) %>%
  left_join(., fat_fit) %>%
  mutate(`r2 (SD)` = ifelse(!model == "Null", sprintf("%.2f (%.2f)", r2, r2_sd), "--"),
         `RMSE (SD)` = ifelse(!model == "Null", sprintf("%.2f (%.2f)", RMSE, RMSE_sd), "--")) %>%
  unite(model, model, mod_sig, sep = "") %>%
  select(model, AIC, `r2 (SD)`, `RMSE (SD)`, term, `b [95% CI]`, SE, t, p) %>%
  arrange(AIC) %>%
  mutate(model = gsub("COMBINED", "Combined", model),
         outcome = "Body fat") %>%
    select(outcome, model, AIC, `r2 (SD)`, `RMSE (SD)`, term, `b [95% CI]`, SE, t, p)

fat_table %>%
  kable(format = "pandoc", digits = 2)

outcome	model	AIC	r2 (SD)	RMSE (SD)	term	b [95% CI]	SE	t	p
Body fat	Null	712.18	–	–	Intercept	0.02 [-0.11, 0.14]	0.06	0.27	.791
Body fat	ROI	712.92	0.02 (0.02)	1.00 (0.05)	Intercept	0.02 [-0.10, 0.15]	0.06	0.38	.706
Body fat	ROI	712.92	0.02 (0.02)	1.00 (0.05)	Reward ROI	-0.14 [-0.39, 0.11]	0.13	-1.12	.265
Body fat	PEV	712.94	0.03 (0.04)	1.00 (0.07)	Intercept	0.00 [-0.12, 0.13]	0.06	0.05	.960
Body fat	PEV	712.94	0.03 (0.04)	1.00 (0.07)	Craving PEV	-0.09 [-0.26, 0.07]	0.08	-1.11	.268
Body fat	Combined	714.49	0.03 (0.03)	1.01 (0.06)	Intercept	0.01 [-0.12, 0.14]	0.07	0.19	.849
Body fat	Combined	714.49	0.03 (0.03)	1.01 (0.06)	Reward ROI	-0.09 [-0.38, 0.19]	0.14	-0.66	.508
Body fat	Combined	714.49	0.03 (0.03)	1.01 (0.06)	Craving PEV	-0.06 [-0.25, 0.13]	0.10	-0.65	.516

EMA enactment

enact_prop_fit = data.frame(model = c("ROI", "PEV", "COMBINED"),
                     r2 = c(mean(best_enact_prop_roi$resample$Rsquared), 
                            mean(best_enact_prop_pev$resample$Rsquared), 
                            mean(enact_prop_combined$resample$Rsquared)),
                     r2_sd = c(sd(best_enact_prop_roi$resample$Rsquared), 
                               sd(best_enact_prop_pev$resample$Rsquared), 
                               sd(enact_prop_combined$resample$Rsquared)),
                     RMSE = c(mean(best_enact_prop_roi$resample$RMSE), 
                              mean(best_enact_prop_pev$resample$RMSE), 
                              mean(enact_prop_combined$resample$RMSE)),
                     RMSE_sd = c(sd(best_enact_prop_roi$resample$RMSE), 
                                 sd(best_enact_prop_pev$resample$RMSE), 
                                 sd(enact_prop_combined$resample$RMSE)))

enact_table = broom::tidy(best_enact_prop_roi$finalModel, conf.int = TRUE) %>% 
  mutate(model = "ROI", 
         mod_sig = ifelse(pf(summary(best_enact_prop_roi$finalModel)$fstatistic[1], 
                             summary(best_enact_prop_roi$finalModel)$fstatistic[2], 
                             summary(best_enact_prop_roi$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", "")) %>%
  bind_rows(broom::tidy(best_enact_prop_pev$finalModel, conf.int = TRUE) %>% 
              mutate(model = "PEV",
                     mod_sig = ifelse(pf(summary(best_enact_prop_pev$finalModel)$fstatistic[1], 
                             summary(best_enact_prop_pev$finalModel)$fstatistic[2], 
                             summary(best_enact_prop_pev$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", ""))) %>%
  bind_rows(broom::tidy(enact_prop_combined$finalModel, conf.int = TRUE) %>% 
              mutate(model = "COMBINED",
                     mod_sig = ifelse(pf(summary(enact_prop_combined$finalModel)$fstatistic[1], 
                             summary(enact_prop_combined$finalModel)$fstatistic[2], 
                             summary(enact_prop_combined$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", ""))) %>%
  bind_rows(broom::tidy(enact_prop_null, conf.int = TRUE) %>% mutate(model = "Null", mod_sig = "")) %>%
  rename("SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%.2f, %.2f]", estimate, conf.low, conf.high),
         p = ifelse(p < .001, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p))),
         term = gsub("\\(Intercept\\)", "Intercept", term),
         term = gsub("_", " ", term),
         term = Hmisc::capitalize(term)) %>%
  left_join(., enact_aic) %>%
  left_join(., enact_prop_fit) %>%
  mutate(`r2 (SD)` = ifelse(!model == "Null", sprintf("%.2f (%.2f)", r2, r2_sd), "--"),
         `RMSE (SD)` = ifelse(!model == "Null", sprintf("%.2f (%.2f)", RMSE, RMSE_sd), "--")) %>%
  unite(model, model, mod_sig, sep = "") %>%
  select(model, AIC, `r2 (SD)`, `RMSE (SD)`, term, `b [95% CI]`, SE, t, p) %>%
  arrange(AIC) %>%
  mutate(model = gsub("COMBINED", "Combined", model),
         outcome = "Enactment") %>%
    select(outcome, model, AIC, `r2 (SD)`, `RMSE (SD)`, term, `b [95% CI]`, SE, t, p)

enact_table %>%
  kable(format = "pandoc", digits = 2)

outcome	model	AIC	r2 (SD)	RMSE (SD)	term	b [95% CI]	SE	t	p
Enactment	PEV*	270.78	0.15 (0.13)	0.96 (0.16)	Intercept	-0.06 [-0.26, 0.13]	0.10	-0.64	.525
Enactment	PEV*	270.78	0.15 (0.13)	0.96 (0.16)	Reward PEV	0.37 [0.14, 0.59]	0.11	3.19	.002
Enactment	Combined*	272.11	0.12 (0.11)	0.97 (0.11)	Intercept	-0.11 [-0.35, 0.12]	0.12	-0.97	.332
Enactment	Combined*	272.11	0.12 (0.11)	0.97 (0.11)	Value ROI	0.13 [-0.19, 0.45]	0.16	0.81	.419
Enactment	Combined*	272.11	0.12 (0.11)	0.97 (0.11)	Reward PEV	0.31 [0.05, 0.58]	0.13	2.33	.022
Enactment	ROI*	275.56	0.11 (0.09)	0.99 (0.08)	Intercept	-0.16 [-0.39, 0.08]	0.12	-1.32	.189
Enactment	ROI*	275.56	0.11 (0.09)	0.99 (0.08)	Value ROI	0.32 [0.04, 0.60]	0.14	2.27	.026
Enactment	Null	278.66	–	–	Intercept	-0.01 [-0.21, 0.19]	0.10	-0.10	.918

combined

bind_rows(bmi_table, fat_table, enact_table) %>%
  kable(format = "pandoc", digits = 2)

outcome	model	AIC	r2 (SD)	RMSE (SD)	term	b [95% CI]	SE	t	p
BMI	Combined*	721.56	0.04 (0.05)	1.00 (0.13)	Intercept	-0.07 [-0.23, 0.08]	0.08	-0.94	.347
BMI	Combined*	721.56	0.04 (0.05)	1.00 (0.13)	Reward ROI	-0.31 [-0.56, -0.05]	0.13	-2.37	.019
BMI	Combined*	721.56	0.04 (0.05)	1.00 (0.13)	Value ROI	0.12 [-0.08, 0.31]	0.10	1.17	.243
BMI	Combined*	721.56	0.04 (0.05)	1.00 (0.13)	Craving regulation PEV	0.14 [-0.05, 0.34]	0.10	1.47	.144
BMI	ROI*	721.74	0.03 (0.05)	1.00 (0.11)	Intercept	-0.03 [-0.17, 0.11]	0.07	-0.43	.667
BMI	ROI*	721.74	0.03 (0.05)	1.00 (0.11)	Reward ROI	-0.31 [-0.56, -0.05]	0.13	-2.35	.020
BMI	ROI*	721.74	0.03 (0.05)	1.00 (0.11)	Value ROI	0.16 [-0.03, 0.35]	0.10	1.61	.108
BMI	PEV	723.45	0.02 (0.03)	0.99 (0.14)	Intercept	-0.05 [-0.19, 0.10]	0.07	-0.67	.507
BMI	PEV	723.45	0.02 (0.03)	0.99 (0.14)	Craving regulation PEV	0.16 [-0.03, 0.34]	0.09	1.64	.102
BMI	Null	724.15	–	–	Intercept	0.01 [-0.11, 0.14]	0.06	0.21	.837
Body fat	Null	712.18	–	–	Intercept	0.02 [-0.11, 0.14]	0.06	0.27	.791
Body fat	ROI	712.92	0.02 (0.02)	1.00 (0.05)	Intercept	0.02 [-0.10, 0.15]	0.06	0.38	.706
Body fat	ROI	712.92	0.02 (0.02)	1.00 (0.05)	Reward ROI	-0.14 [-0.39, 0.11]	0.13	-1.12	.265
Body fat	PEV	712.94	0.03 (0.04)	1.00 (0.07)	Intercept	0.00 [-0.12, 0.13]	0.06	0.05	.960
Body fat	PEV	712.94	0.03 (0.04)	1.00 (0.07)	Craving PEV	-0.09 [-0.26, 0.07]	0.08	-1.11	.268
Body fat	Combined	714.49	0.03 (0.03)	1.01 (0.06)	Intercept	0.01 [-0.12, 0.14]	0.07	0.19	.849
Body fat	Combined	714.49	0.03 (0.03)	1.01 (0.06)	Reward ROI	-0.09 [-0.38, 0.19]	0.14	-0.66	.508
Body fat	Combined	714.49	0.03 (0.03)	1.01 (0.06)	Craving PEV	-0.06 [-0.25, 0.13]	0.10	-0.65	.516
Enactment	PEV*	270.78	0.15 (0.13)	0.96 (0.16)	Intercept	-0.06 [-0.26, 0.13]	0.10	-0.64	.525
Enactment	PEV*	270.78	0.15 (0.13)	0.96 (0.16)	Reward PEV	0.37 [0.14, 0.59]	0.11	3.19	.002
Enactment	Combined*	272.11	0.12 (0.11)	0.97 (0.11)	Intercept	-0.11 [-0.35, 0.12]	0.12	-0.97	.332
Enactment	Combined*	272.11	0.12 (0.11)	0.97 (0.11)	Value ROI	0.13 [-0.19, 0.45]	0.16	0.81	.419
Enactment	Combined*	272.11	0.12 (0.11)	0.97 (0.11)	Reward PEV	0.31 [0.05, 0.58]	0.13	2.33	.022
Enactment	ROI*	275.56	0.11 (0.09)	0.99 (0.08)	Intercept	-0.16 [-0.39, 0.08]	0.12	-1.32	.189
Enactment	ROI*	275.56	0.11 (0.09)	0.99 (0.08)	Value ROI	0.32 [0.04, 0.60]	0.14	2.27	.026
Enactment	Null	278.66	–	–	Intercept	-0.01 [-0.21, 0.19]	0.10	-0.10	.918

VIF

# BMI
car::vif(bmi_combined$finalModel)

##             reward_ROI              value_ROI craving_regulation_PEV 
##               1.101738               1.185011               1.085275

# %fat
car::vif(fat_combined$finalModel)

##  reward_ROI craving_PEV 
##    1.303633    1.303633

# enactment
car::vif(enact_prop_combined$finalModel)

##  value_ROI reward_PEV 
##   1.348297   1.348297

Visualization & Cross-validated Model Comparison

Dani Cosme & Rich Lopez

2020-09-17

load packages

define palettes

source data

standardize brain data

summarize outliers

merge data

tidy data for plotting

tidy data for modeling

visualize

DVs and dietary restraint

combined

ROIs

means

correlations

PEVs

means

correlations

correlations among ROIs and PEVs

PEVs v ROIs

ROIs

PEVs

descriptives

means and sds

correlations

t-tests

run cross-validated models

specify model variables

BMI

ROI models

PEV models

combined

summarize the overall best fitting

% fat

ROI models

PEV models

combined

summarize the overall best fitting

EMA enactment proportion

ROI models

PEV models

combined

summarize the overall best fitting

tables

bmi

% fat

EMA enactment

combined

VIF